home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / gravsim.com / GRAVSIM.C < prev   
Encoding:
C/C++ Source or Header  |  1989-09-13  |  14.2 KB  |  376 lines
  1. /************************************************************************
  2.           2-Dimensional Planetary Simulation for Turbo C
  3.                by C. Daniel Hassell
  4.  
  5. Adapted from "Force-Based Simulations" by Todd King in Dr. Dobb's Journal,
  6. September 1989.  The original code was in C++, and operated in character
  7. mode with no run-time controls.  (No criticism of Mr. King intended; he only
  8. wanted to demo the power of C++ in a short program.  I don't own C++ and
  9. don't want to, but I liked his simulation concept.)
  10.  
  11. The program simulates the movements of a set of bodies affecting each other
  12. by mutual gravitational pull.  The simulation is quite realistic, except for
  13. its 2-dimentional nature.  The masses are in Kg and the velocities are in
  14. Km/Hr.  For example, the entire list of 9 planets and their moons could be
  15. included and shown orbiting the Sun.
  16.  
  17. The view space is defined by extent_x and extent_y.  Zooming out doubles the
  18. values of these variables, while zooming in halves them.  A circle is drawn
  19. for each planet within the view space with a radius related to the cube root
  20. of the mass of the planet.  Panning and zooming in or out is allowed using
  21. the PC cursor movement keys.  Panning in any direction occurs in 1/4-screen
  22. increments.
  23.  
  24. Each loop through the simulation equals one hour.  The time is shown in the
  25. upper left corner of the screen as days:hours.  The speed of the simulation
  26. is dependent on the speed of the system; it can be slowed down by adding a
  27. delay() function somewhere but not easily cranked up faster.  On my old
  28. PC/AT with a math chip, it takes a little over a second to tick off one
  29. simulated day when there are two bodies in the system.  Each added planet
  30. slows it down noticably, since the number of FLOPS required is a multiple
  31. of n * ( n - 1 ), where n = the number of bodies.
  32.  
  33. The initial set of planets must now be entered as constants in the main()
  34. function, but I hope some bright person designs a better user interface.
  35. (Turbo C's easy environment can make one lazy that way.)  A new "planet X",
  36. defined in the big_bang() function, can be dropped in by pressing 'x'.
  37.  
  38. Collisions between planets result in the absorption of the two into one.
  39. The resulting body has a mass and velocity vector equal to the sum of those
  40. of the two separate bodies.  The total number of planets existing at one
  41. time is limited to POOL_LIMIT.  On a color display, the planets will
  42. leave a trail of blue so that it is easy to see where they have been, but
  43. the trails are lost when the pan or zoom functions are used.  On a monochrome
  44. monitor, no trail is left.
  45.  
  46. The simulation is now set up as it was published by Todd King, showing the
  47. Moon orbiting the Earth.  A "planet X" can be introduced at any time to
  48. demonstrate the disturbing effect a passing body could have on the orbits
  49. of the Earth and the Moon.  Any hypothetical orbital system can be written
  50. in, however.  An interesting arrangement is to set four identical planets
  51. at the corners of a perfect square in orbit around their common center of
  52. gravity.  No matter how exact the initial placement, this system is highly
  53. unstable and will become chaotic over time.
  54.  
  55. Comments are welcome at [72017,2466].
  56.  
  57. C. Daniel Hassell,  September 13, 1989.
  58.  
  59. ************************************************************************/
  60. #include <dos.h>
  61. #include <stdio.h>
  62. #include <stdlib.h>
  63. #include <conio.h>
  64. #include <math.h>
  65. #include <time.h>
  66. #include <graphics.h>
  67.  
  68. #define G            -6.67e-11     /* the gravitational constant */
  69. #define SECS_PER_TIC  3600L        /* 1 Hour, originally = 86400 or 1 day */
  70. #define ESC           27
  71. #define POOL_LIMIT 20           /* limit on the number of bodies */
  72. #define RADIUS_FUDGE  50           /* factor used to set size of circles */
  73.  
  74. typedef struct {
  75.   double x;
  76.   double y;
  77. } VECTOR_2D;                /* vector typedef for location and velocity */
  78.  
  79. typedef struct {
  80.   VECTOR_2D world;          /* location in Km */
  81.   VECTOR_2D velocity;       /* velocity in Km/Hr */
  82.   double mass;              /* mass in Kg */
  83.   double gmass;             /* G * mass */
  84.   int radius;               /* Radius of circle for this body */
  85. } BODY;
  86.  
  87. /* Global declarations */
  88. unsigned body_cnt;              /* BODY count (pardon the expression) */
  89. BODY   *body_pool[POOL_LIMIT];  /* array of pointers to BODYs */
  90. int    GraphDriver;        /* the graphics device driver        */
  91. int    GraphMode;        /* the graphics mode value        */
  92. int    MaxX, MaxY;        /* the maximum resolution of the screen */
  93. int    MaxColor;        /* the maximum color # available    */
  94. int    ErrorCode;        /* reports any graphics errors        */
  95. int    height;                  /* height of default text               */
  96. double extent_x = 15e8;
  97. double extent_y = 15e8;         /* initial size of view space in Km */
  98. int    trail_color;             /* color of trail behind planets */
  99.  
  100. /* Function prototypes */
  101. void initialize(void);
  102. void drawborder(void);
  103. int gprintf(int xloc, int yloc, char *fmt, ... );
  104. void set_mass( BODY *b, double m );
  105. void set_velocity( BODY *b, double x, double y );
  106. void set_position( BODY *b, double x, double y );
  107. void apply_force( BODY *b, BODY *from );
  108. void update( BODY *b );
  109. void print_tick( void );
  110. void convert( VECTOR_2D *world, VECTOR_2D *screen );
  111. void collision( BODY *b, BODY *from );
  112. void make_planet( double mass, double px, double py, double vx, double vy );
  113. char move_view( void );
  114. void big_bang( void );
  115.  
  116.  
  117. void initialize(void)
  118. { /*  Initializes the graphics system and reports any errors */
  119.   GraphDriver = DETECT;         /* Request auto-detection    */
  120.   initgraph( &GraphDriver, &GraphMode, "" );
  121.   ErrorCode = graphresult();        /* Read result of initialization*/
  122.   if( ErrorCode != grOk ){        /* Error occured during init    */
  123.     printf(" Graphics System Error: %s\n", grapherrormsg( ErrorCode ) );
  124.     exit( 1 );
  125.   }
  126.   MaxColor = getmaxcolor();        /* Read maximum number of colors*/
  127.   MaxX = getmaxx();
  128.   MaxY = getmaxy();            /* Read size of screen        */
  129.   return;
  130. }
  131.  
  132. void drawborder(void)
  133. { /*  Draw a solid line border and establish text settings. */
  134.   cleardevice();            /* Clear graphics screen    */
  135.   height = textheight( "H" );           /* Get basic text height        */
  136.   setcolor( MaxColor );
  137.   setlinestyle( SOLID_LINE, 0, NORM_WIDTH );
  138.   rectangle( 0, height+4, MaxX, MaxY-(height+4) );    /* Draw border */
  139.   settextjustify( CENTER_TEXT, TOP_TEXT );
  140.   gprintf( MaxX/2, 1, "Planetary Simulation" );
  141.   settextjustify( CENTER_TEXT, BOTTOM_TEXT );
  142.   gprintf( MaxX/2, MaxY, "Press ESC to quit, PgUp or PgDn to zoom, X for new planet, or %c%c%c%c to pan.", 24, 25, 26, 27 );
  143.   settextjustify( LEFT_TEXT, TOP_TEXT );
  144.   if ( MaxColor > 2 )       /* set color of trail left by planets */
  145.     trail_color = BLUE;     /* blue trail if colors available */
  146.   else
  147.     trail_color = BLACK;    /* ... else no trail */
  148.   return;
  149. }
  150.  
  151. int gprintf( int xloc, int yloc, char *fmt, ... )
  152. { /* Used like PRINTF except the output is sent to the
  153.      screen in graphics mode at the specified co-ordinate.            */
  154.   va_list  argptr;            /* Argument list pointer    */
  155.   char str[140];            /* Buffer to build sting into    */
  156.   int cnt;                /* Result of SPRINTF for return */
  157.  
  158.   va_start( argptr, format );        /* initialize va_ functions    */
  159.   cnt = vsprintf( str, fmt, argptr );    /* prints string to buffer    */
  160.   outtextxy( xloc, yloc, str );            /* Send string in graphics mode */
  161.   va_end( argptr );            /* Close va_ functions        */
  162.   return( cnt );            /* Return the conversion count    */
  163. }
  164.  
  165. void set_mass( BODY *b, double m )
  166. {  /* radius determined by cube root of mass */
  167.    b->mass = m;
  168.    b->gmass = G * m;
  169.    b->radius = max( 1, (int)( pow(b->mass, (1.0/3.0))/extent_x*RADIUS_FUDGE ));
  170.    return;
  171. }
  172.  
  173. void set_velocity( BODY *b, double x, double y )
  174. {
  175.    b->velocity.x = x;
  176.    b->velocity.y = y;
  177.    return;
  178. }
  179.  
  180. void set_position( BODY *b, double x, double y )
  181. {
  182.    b->world.x = x;
  183.    b->world.y = y;
  184.    return;
  185. }
  186.  
  187. void collision( BODY *b, BODY *from )
  188. { /* Simulates a collision by absorbing from into b */
  189.   double x, y, new_mass;
  190.   int i;
  191.   VECTOR_2D s;
  192.  
  193.   new_mass = b->mass + from->mass;      /* calculate combined mass & speed */
  194.   x = ( b->velocity.x * b->mass + from->velocity.x * from->mass ) / new_mass;
  195.   y = ( b->velocity.y * b->mass + from->velocity.y * from->mass ) / new_mass;
  196.   convert( &from->world, &s );
  197.   setcolor( trail_color );
  198.   circle( (int)s.x, (int)s.y, from->radius );  /* removes from from screen */
  199.   convert( &b->world, &s );
  200.   circle( (int)s.x, (int)s.y, b->radius );  /* removes b from screen */
  201.   set_mass( b, new_mass );
  202.   set_velocity( b, x, y );
  203.   setcolor( MaxColor );
  204.   circle( (int)s.x, (int)s.y, b->radius );  /* put new body on screen */
  205.   for ( i = 0; i < (body_cnt-1); i++ ) {
  206.     if ( body_pool[ i ] == from ) {
  207.       body_pool[ i ] = body_pool[ i+1 ];
  208.       body_pool[ i+1 ] = from;       /* shift other elements of array down */
  209.     }                                /*  ... to fill gap left by absorbed body */
  210.   }                                  /* pointer to be eliminated now at top */
  211.   free( body_pool[ --body_cnt ] );   /* free top pointer in array */
  212.   body_pool[ body_cnt ] = NULL;      /* set pointer to NULL to prevent confusion */
  213.   return;
  214. }
  215.  
  216. void apply_force( BODY *b, BODY *from )
  217. { /* Adjusts velocity vector of BODY b given the force of gravity from FROM */
  218.   VECTOR_2D d;
  219.   double rs, r, v, scrn_dist;
  220.  
  221.   if ( (b == from) || (b == NULL) || (from == NULL) )
  222.     return;
  223.   d.x = b->world.x - from->world.x;
  224.   d.y = b->world.y - from->world.y;     /* d vector now the distance between bodies */
  225.   rs = ( d.x * d.x ) + ( d.y * d.y );   /* finds straight line distance */
  226.   if ( rs != 0.0 ) {                    /* ... by calculating hypotenuse */
  227.     r = sqrt( rs );                     /* now r = straight line distance */
  228.     scrn_dist = MaxX * r / extent_x;    /* scrn_dist is in pixels */
  229.     if ( scrn_dist < ( b->radius + from->radius )) {
  230.       collision( b, from );             /* if the circles touch, call collision() */
  231.     } else {
  232.       v = ( from->gmass / rs ) * SECS_PER_TIC;  /* v = the change in velocity */
  233.       b->velocity.x += v * d.x / r;     /* add to x and y elements of vector */
  234.       b->velocity.y += v * d.y / r;
  235.     }
  236.   }
  237.   return;
  238. }
  239.  
  240. void convert( VECTOR_2D *world, VECTOR_2D *screen )
  241. {  /* Converts world coordinates to screen coordinates */
  242.    screen->x = MaxX * ( world->x / extent_x );
  243.    screen->y = MaxY * ( world->y / extent_y );
  244.    return;
  245. }
  246.  
  247. void update( BODY *b )
  248. { /* Updates the position coordinates of BODY b and moves circle on screen */
  249.    VECTOR_2D s;
  250.  
  251.    convert( &b->world, &s );
  252.    if (( s.x < MaxX ) && ( s.x >= 0 ) && ( s.y < MaxY ) && ( s.y >= 0 )) {
  253.      setcolor( trail_color );    /* erase old circle */
  254.      circle( (int)s.x, (int)s.y, b->radius );
  255.    }
  256.    b->world.x += b->velocity.x * SECS_PER_TIC;      /* update world coord. */
  257.    b->world.y += b->velocity.y * SECS_PER_TIC;
  258.    convert( &b->world, &s );
  259.    if (( s.x < MaxX ) && ( s.x >= 0 ) && ( s.y < MaxY ) && ( s.y >= 0 )) {
  260.      setcolor( MaxColor );        /* draw new circle */
  261.      circle( (int)s.x, (int)s.y, b->radius );
  262.    }
  263.    return;
  264. }
  265.  
  266. void print_tick( void )
  267. {  /* Prints day:hour in upper left corner */
  268.    static unsigned int tick = 0;
  269.  
  270.    setcolor( BLACK );
  271.    gprintf(1, 1, "%4d:%2u ", (tick / 24), (tick % 24) ); /* erase old one */
  272.    tick++;
  273.    setcolor( MaxColor );
  274.    gprintf(1, 1, "%4d:%2u ", (tick / 24), (tick % 24) ); /* write new one */
  275.    return;
  276. }
  277.  
  278. void make_planet( double mass, double px, double py, double vx, double vy )
  279. {  /* creates a new BODY and assigns the next available pointer to it */
  280.    BODY *newplanet;
  281.  
  282.    if ( body_cnt >= POOL_LIMIT )
  283.     return;
  284.    newplanet = calloc( 1, sizeof( BODY ) );
  285.    set_mass( newplanet, mass );
  286.    set_position( newplanet, px, py );
  287.    set_velocity( newplanet, vx, vy );
  288.    body_pool[ body_cnt++ ] = newplanet;
  289.    return;
  290. }
  291.  
  292. char move_view( void )
  293. {  /* interprets keystrokes and moves view space accordingly */
  294.    char r = 0;
  295.    int i;
  296.  
  297.    if ( !kbhit() )
  298.      return( r );
  299.    switch ( r = getch() ) {
  300.      case 'X' :        /* If X or x is pressed then release planet x */
  301.      case 'x' :
  302.        make_planet( 5e23, extent_x, 0, -1600, 2000);
  303.        break;
  304.      case ESC :        /* ESC to quit */
  305.        break;
  306.      case 0:
  307.        switch ( getch() ) {
  308.        case 73:    /* PgUp */
  309.      extent_x *= 2;
  310.      extent_y *= 2;    /* zoom out */
  311.      for ( i = 0; i < body_cnt; i++ ) {
  312.        set_mass( body_pool[i], body_pool[i]->mass );  /* reset radii */
  313.        body_pool[i]->world.x += extent_x/4;    /* keep view centered */
  314.        body_pool[i]->world.y += extent_y/4;
  315.      }
  316.      break;
  317.        case 81:    /* PgDn */
  318.      extent_x /= 2;
  319.      extent_y /= 2;    /* zoom in */
  320.      for ( i = 0; i < body_cnt; i++ ) {
  321.        set_mass( body_pool[i], body_pool[i]->mass );  /* reset radii */
  322.        body_pool[i]->world.x -= extent_x/2;   /* keep view centered */
  323.        body_pool[i]->world.y -= extent_y/2;
  324.      }
  325.      break;
  326.        case 72:    /* UP arrow */
  327.      for ( i = 0; i < body_cnt; i++ )
  328.        body_pool[i]->world.y += extent_y/4;
  329.      break;
  330.        case 75:    /* LEFT arrow */
  331.      for ( i = 0; i < body_cnt; i++ )
  332.        body_pool[i]->world.x += extent_x/4;
  333.      break;
  334.        case 77:    /* RIGHT arrow */
  335.      for ( i = 0; i < body_cnt; i++ )
  336.        body_pool[i]->world.x -= extent_x/4;
  337.      break;
  338.        case 80:    /* DOWN arrow */
  339.      for ( i = 0; i < body_cnt; i++ )
  340.        body_pool[i]->world.y -= extent_y/4;
  341.      break;
  342.      }  /* end of extended char switch */
  343.      drawborder();    /* clear screen and redraw border */
  344.    }  /* end of simple char switch */
  345.    return( r );
  346. }
  347.  
  348. void big_bang( void )
  349. {  /* Main loop of simulation */
  350.    int i, j;
  351.  
  352.    do {
  353.      print_tick();
  354.      for ( i = 0; i < body_cnt; i++ ) {
  355.        for ( j = 0; j < body_cnt; j++ )   /* apply to each pair of bodies */
  356.      apply_force( body_pool[ i ], body_pool[ j ] );
  357.      }
  358.      for ( i = 0; i < body_cnt; i++ )   /* update locations */
  359.        update( body_pool[ i ] );
  360.    } while ( move_view() != ESC );
  361.    return;                              /* quit if ESC was pressed */
  362. }
  363.  
  364. void main( void )
  365. {
  366.   initialize();
  367.   drawborder();
  368.  
  369.   make_planet( 5.9e24,  7.5e8,  7.5e8,  12.5,     0 );
  370.   make_planet( 7.4e22,  7.5e8, 11.3e8, -1020,     0 );
  371.   big_bang();
  372.  
  373.   closegraph();
  374. }
  375.  
  376.